This document is for preliminary exploration of the data from the Subtidal Ecology lab at Pontificia Universidad. See data-collection.html for details on the collection.
Start with the simplest. Holdfast diameter will be the primary metric for tracking plant size and will be used to estimate everything else, including biomass and demographic rates.
morph.df <- list(
IFOP_morph=read_csv(glue("{data.dir}IFOP_Kelp_Morphology.csv"),
na=c("NaN", "")) %>%
select(c(1:9,11:13,15:19)) %>%
filter(Group=="Macroalga" &
Species=="Lessonia trabeculata" &
Abundance > 0) %>%
mutate(Date=ymd(paste(Year, Month, Day, sep="-")),
Site=factor(Site),
Management=factor(Management),
Frondosity=factor(Frondosity),
source="IFOP_morph") %>%
rename(LengthTot=TL),
IFOP_allom=read_csv(glue("{data.dir}IFOP_allometry.csv")) %>%
select(-16) %>%
mutate(Date=ymd(paste(Year, Month, Day, sep="-")),
Site=factor(Site),
Management=factor(Management),
Weight=1e3*Weight,
source="IFOP_allom") %>%
rename(LengthTot=TL,
WeightTot=Weight),
NERC_bl=read_csv(glue("{data.dir}Quads_baseline.csv")) %>%
select(-c(1,3:4,10:12)) %>%
filter(Species=="Lessonia trabeculata") %>%
mutate(Site=factor(Site),
Zone=factor(Zone),
Management=factor(Management),
Upwelling=factor(Upwelling),
Frondosity=factor(Frondosity),
source="NERC_bl") %>%
rename(Holdfast=Holdfast.Diameter,
LengthTot=Total.Length),
NERC_ph=read_csv(glue("{data.dir}Patches_Morphology.csv")) %>%
mutate(Site=factor(Site),
Zone=factor(Zone),
Management=factor(Management),
Upwelling=factor(Upwelling),
source="NERC_ph") %>%
rename(Holdfast=Holdfast.Diameter,
LengthTot=Total.Length),
Cata_ph=read_csv(glue("{data.dir}Morfo_DW_C.csv")) %>%
select(-c(1,11,20:30)) %>%
mutate(Site=factor(Site),
Zone=factor(Zone),
Management=factor(Management)) %>%
rename(Holdfast=Diam_Holdfast,
Height=High_Holdfast,
LengthTot=LT,
WeightTot=TW_sum,
Weight=Weight_Str,
Length=LT_stipes,
Dichot=Dichot_stipes,
Diam=Diam_stipes) %>%
pivot_wider(names_from="Structure", values_from=c(8:17)) %>%
select(!where(~all(is.na(.x)))) %>%
select(-starts_with("Structure_")) %>%
mutate(source="Cata_ph") %>%
rename(Stipes=Stipes_Stipes)
) %>%
do.call('bind_rows', .)
glimpse(morph.df)
## Rows: 4,563
## Columns: 49
## $ Serie <dbl> 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 1…
## $ Date <date> 2018-05-16, 2018-05-16, 2018-05-16, 2018-05-16, 20…
## $ Day <dbl> 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,…
## $ Month <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ Year <dbl> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 201…
## $ Site <fct> Quintay A, Quintay A, Quintay A, Quintay A, Quintay…
## $ Management <fct> TURF, TURF, TURF, TURF, TURF, TURF, TURF, TURF, TUR…
## $ Transect <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ Quadrats <dbl> 1, 2, 2, 3, 4, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 8, …
## $ Depth <dbl> 20.0, 18.8, 18.8, 19.0, 18.7, 18.8, 18.8, NA, NA, N…
## $ Group <chr> "Macroalga", "Macroalga", "Macroalga", "Macroalga",…
## $ Species <chr> "Lessonia trabeculata", "Lessonia trabeculata", "Le…
## $ Abundance <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ Holdfast <dbl> 33, 16, 25, 23, 23, 37, 38, 25, 10, 7, 7, 22, 12, 2…
## $ Stipes <dbl> 16, 4, 3, 6, 7, 7, 7, 7, 5, 3, 14, 12, 3, 15, 10, 1…
## $ LengthTot <dbl> 128, 179, 133, 160, 118, 171, 213, 60, 50, 100, 100…
## $ Frondosity <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ source <chr> "IFOP_morph", "IFOP_morph", "IFOP_morph", "IFOP_mor…
## $ Extraction <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ WeightTot <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Holdfast_perimeter <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Reproductive <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Zone <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Upwelling <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Station <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Quad <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ ...1 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Survey.Number <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MAE <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Latitude <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Longitude <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Surveyor.Name <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Patch <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Replicate <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ ID <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Holdfast <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Weight_Blades <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Holdfast <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FW_Blades <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Holdfast <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DW_Blades <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Height_Holdfast <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ First_bif_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Diam_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Length_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Dichot_Stipes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
summary(morph.df)
## Serie Date Day Month
## Min. : 1.0 Min. :2018-04-20 Min. : 2.00 Min. : 1.00
## 1st Qu.: 993.8 1st Qu.:2019-11-10 1st Qu.: 7.00 1st Qu.: 4.00
## Median :2618.5 Median :2020-10-04 Median :16.00 Median : 8.00
## Mean :2257.7 Mean :2020-06-12 Mean :14.91 Mean : 7.22
## 3rd Qu.:3407.2 3rd Qu.:2021-04-19 3rd Qu.:24.00 3rd Qu.:11.00
## Max. :4240.0 Max. :2022-05-17 Max. :30.00 Max. :12.00
## NA's :2131 NA's :87 NA's :2131 NA's :2131
## Year Site Management
## Min. :2018 Isla Grande de Atacama : 503 OA :1037
## 1st Qu.:2018 La Ballena : 382 Reserve: 154
## Median :2020 Chanaral de Aceituno OA : 370 TURF :3027
## Mean :2020 Puerto Viejo : 311 AMERB : 345
## 3rd Qu.:2021 Cobquecura : 287
## Max. :2022 Chanaral de Aceituno TURF: 280
## NA's :2131 (Other) :2430
## Transect Quadrats Depth Group
## Min. : 1.000 Min. : 1.000 Min. : 3.30 Length:4563
## 1st Qu.: 2.000 1st Qu.: 3.000 1st Qu.: 9.10 Class :character
## Median : 3.000 Median : 4.000 Median :11.00 Mode :character
## Mean : 4.167 Mean : 4.537 Mean :11.21
## 3rd Qu.: 6.000 3rd Qu.: 6.000 3rd Qu.:12.80
## Max. :10.000 Max. :10.000 Max. :24.60
## NA's :1723 NA's :2451 NA's :1806
## Species Abundance Holdfast Stipes
## Length:4563 Min. :0.0000 Min. : 0.00 Min. : 0.000
## Class :character 1st Qu.:1.0000 1st Qu.: 7.00 1st Qu.: 2.000
## Mode :character Median :1.0000 Median :19.00 Median : 3.000
## Mean :0.9996 Mean :19.28 Mean : 5.091
## 3rd Qu.:1.0000 3rd Qu.:29.00 3rd Qu.: 6.667
## Max. :1.0000 Max. :95.00 Max. :79.000
## NA's :1723 NA's :1 NA's :86
## LengthTot Frondosity source Extraction
## Min. : 0.0 1 : 56 Length:4563 Length:4563
## 1st Qu.: 50.0 2 : 358 Class :character Class :character
## Median :131.0 3 :1621 Mode :character Mode :character
## Mean :126.6 0 : 1
## 3rd Qu.:193.0 NA's:2527
## Max. :450.0
## NA's :18
## WeightTot Holdfast_perimeter Reproductive Zone
## Min. : 1.0 Min. : 0.80 Min. :0.000 Central : 606
## 1st Qu.: 332.5 1st Qu.: 9.00 1st Qu.:0.000 North : 814
## Median : 1380.0 Median : 18.00 Median :0.000 North-Central: 650
## Mean : 2637.0 Mean : 20.23 Mean :0.424 Center : 28
## 3rd Qu.: 3936.5 3rd Qu.: 27.50 3rd Qu.:1.000 North-Center : 33
## Max. :20000.0 Max. :161.00 Max. :1.000 NA's :2432
## NA's :4156 NA's :4376 NA's :4280
## Upwelling Station Quad ...1
## nonupwelling: 503 Min. : 0.00 Min. : 1.000 Min. : 1.0
## upwelling : 225 1st Qu.: 20.00 1st Qu.: 3.000 1st Qu.: 329.8
## N : 891 Median : 40.00 Median : 5.000 Median : 658.5
## Y : 425 Mean : 41.05 Mean : 4.919 Mean : 658.5
## NA's :2519 3rd Qu.: 60.00 3rd Qu.: 7.000 3rd Qu.: 987.2
## Max. :100.00 Max. :10.000 Max. :1316.0
## NA's :3835 NA's :3835 NA's :3247
## Survey.Number MAE Latitude Longitude
## Min. :0.000 Min. : 0.00 Min. :-32.31 Min. :-71.53
## 1st Qu.:1.000 1st Qu.:12.00 1st Qu.:-32.21 1st Qu.:-71.50
## Median :3.000 Median :24.00 Median :-29.05 Median :-71.48
## Mean :2.337 Mean :19.22 Mean :-29.41 Mean :-71.30
## 3rd Qu.:3.000 3rd Qu.:24.00 3rd Qu.:-27.34 3rd Qu.:-70.98
## Max. :4.000 Max. :30.00 Max. :-27.25 Max. :-70.96
## NA's :3247 NA's :3247 NA's :3257 NA's :3257
## Surveyor.Name Patch Replicate ID
## Length:4563 Min. :1.000 Min. : 1.000 Min. : 1.000
## Class :character 1st Qu.:1.000 1st Qu.: 4.000 1st Qu.: 4.000
## Mode :character Median :2.000 Median : 7.000 Median : 8.000
## Mean :2.438 Mean : 8.362 Mean : 8.138
## 3rd Qu.:3.000 3rd Qu.:12.000 3rd Qu.:11.000
## Max. :4.000 Max. :26.000 Max. :22.000
## NA's :3247 NA's :3252 NA's :4476
## Weight_Holdfast Weight_Stipes Weight_Blades FW_Holdfast
## Min. : 1.0 Min. : 0.0 Min. : 14 Min. : 0.179
## 1st Qu.: 46.0 1st Qu.: 53.0 1st Qu.: 187 1st Qu.: 1.103
## Median : 178.0 Median : 377.0 Median : 471 Median : 1.992
## Mean : 390.2 Mean : 831.4 Mean :1029 Mean : 2.595
## 3rd Qu.: 575.5 3rd Qu.:1152.0 3rd Qu.:1088 3rd Qu.: 3.245
## Max. :3200.0 Max. :8615.0 Max. :6885 Max. :13.533
## NA's :4476 NA's :4476 NA's :4476 NA's :4476
## FW_Stipes FW_Blades DW_Holdfast DW_Stipes
## Min. :0.191 Min. :0.023 Min. :0.040 Min. :0.045
## 1st Qu.:1.050 1st Qu.:0.526 1st Qu.:0.289 1st Qu.:0.224
## Median :1.891 Median :0.666 Median :0.492 Median :0.485
## Mean :2.311 Mean :0.662 Mean :0.784 Mean :0.673
## 3rd Qu.:3.055 3rd Qu.:0.821 3rd Qu.:0.968 3rd Qu.:0.792
## Max. :7.161 Max. :1.529 Max. :4.041 Max. :3.332
## NA's :4477 NA's :4476 NA's :4476 NA's :4477
## DW_Blades Height_Holdfast First_bif_Stipes Diam_Stipes
## Min. :0.030 Min. : 1.00 Min. : 0.50 Min. : 1.00
## 1st Qu.:0.107 1st Qu.: 5.00 1st Qu.: 2.65 1st Qu.: 12.17
## Median :0.132 Median : 11.00 Median : 8.10 Median : 26.60
## Mean :0.140 Mean : 12.74 Mean :15.93 Mean : 27.97
## 3rd Qu.:0.177 3rd Qu.: 17.00 3rd Qu.:26.17 3rd Qu.: 37.20
## Max. :0.446 Max. :101.00 Max. :79.00 Max. :296.00
## NA's :4476 NA's :4476 NA's :4476 NA's :4478
## Length_Stipes Dichot_Stipes
## Min. : 0.933 Min. : 0.933
## 1st Qu.: 19.175 1st Qu.: 19.175
## Median : 62.500 Median : 62.500
## Mean : 63.080 Mean : 63.080
## 3rd Qu.: 97.500 3rd Qu.: 97.500
## Max. :187.333 Max. :187.333
## NA's :4479 NA's :4479
ggplot(morph.df, aes(Holdfast)) + geom_histogram() +
facet_wrap(~source, scales="free_y")
ggplot(morph.df, aes(Holdfast, colour=source)) + geom_density()
ggplot(morph.df, aes(LengthTot)) + geom_histogram() +
facet_wrap(~source, scales="free_y")
ggplot(morph.df, aes(LengthTot, colour=source)) + geom_density()
ggplot(morph.df, aes(WeightTot)) + geom_histogram() +
facet_wrap(~source, scales="free_y")
ggplot(morph.df, aes(WeightTot, colour=source)) + geom_density()
ggplot(morph.df, aes(Frondosity)) + geom_bar() +
facet_wrap(~source, scales="free_y")
ggplot(morph.df, aes(factor(Reproductive))) + geom_bar() +
facet_wrap(~source, scales="free_y")
ggplot(morph.df, aes(Holdfast, Stipes, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(Holdfast, LengthTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site),
formula=y~x+I(x^2)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source),
formula=y~x+I(x^2)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(Holdfast, WeightTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site),
formula=y~x+I(x^2)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source),
formula=y~x+I(x^2)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(LengthTot, WeightTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous(trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(Holdfast, Weight_Holdfast, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(Holdfast, Weight_Stipes, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
ggplot(morph.df, aes(Holdfast, Weight_Blades, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous(trans="log1p")
# Larger holdfast = less blade-y by weight, more holdfast-y & stipe-y
ggplot(morph.df, aes(Holdfast, Weight_Holdfast/WeightTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous("Holdfast weight / Total", limits=c(0,1))
ggplot(morph.df, aes(Holdfast, Weight_Stipes/WeightTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous("Stipe weight / Total", limits=c(0,1))
ggplot(morph.df, aes(Holdfast, Weight_Blades/WeightTot, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.5, linetype=3, aes(group=Site)) +
stat_smooth(method="lm", se=F, size=0.75, aes(group=source)) +
scale_x_continuous("Holdfast diameter", trans="log1p") +
scale_y_continuous("Blade weight / Total", limits=c(0,1))
ggplot(morph.df, aes(Holdfast, Height_Holdfast, colour=source)) +
geom_point(alpha=0.5, shape=1) +
stat_smooth(method="lm", se=F, size=0.25, aes(group=Site)) +
scale_x_log10() + scale_y_log10()
ggplot(morph.df, aes(Holdfast, colour=factor(Frondosity))) +
geom_density() + scale_colour_viridis_d() +
scale_x_log10()
# Lots of NAs -- Cata's thesis should have this too
ggplot(morph.df, aes(Holdfast, colour=factor(Reproductive))) +
geom_density() +
scale_x_log10()
Growth_Surv <- read_csv(glue("{data.dir}Growth&Survival.csv")) %>%
mutate(Site=factor(Site),
Management=factor(Management),
Plant=factor(Plant),
alive=as.numeric(Holdfast > 0),
Date_Harvest=Date_harvest_dmy %>%
str_replace("Dic", "01-12") %>%
str_replace("Nov", "01-11") %>%
str_replace("Ene", "01-01") %>%
str_replace("Feb", "01-02") %>%
dmy(),
Date_Obs=Date_monitoring %>%
str_replace("Abr", "01-04") %>%
str_replace("Oct", "01-10") %>%
str_replace("May", "01-05") %>%
dmy()) %>%
mutate(days_since_harvest=as.numeric(Date_Obs - Date_Harvest))
survey.i <- Growth_Surv %>%
group_by(Site, Management, Patch, Monitoring) %>%
summarise(Date_Harvest=min(Date_Harvest),
Date_Obs=min(Date_Obs)) %>%
arrange(Management, Site, Patch, Monitoring)
survey.i <- bind_rows(
survey.i,
survey.i %>% slice_head(n=1) %>%
mutate(Monitoring=0,
Date_Obs=Date_Harvest)
) %>%
arrange(Management, Site, Patch, Monitoring) %>%
mutate(days_since_harvest=as.numeric(Date_Obs - Date_Harvest),
days_since_prevObs=as.numeric(Date_Obs - lag(Date_Obs)),
days_since_prevObs=replace_na(days_since_prevObs, 0),
days_to_nextObs=as.numeric(lead(Date_Obs) - Date_Obs))
We can link probability of survival to holdfast size x Management. With the magic of Hierarchical Bayesian modelling, we can annualize the estimated survival rates if we assume that the survival probability is constant throughout the year.
surv.df <- Growth_Surv %>%
select(Site, Management, Patch, Monitoring,
Plant, Holdfast, Stipes, Total_Length) %>%
left_join(., survey.i) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
mutate(prevSurv=Monitoring-1) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
group_by(Site, Patch, Plant) %>%
mutate(doesSurvive=lead(Holdfast>0),
maxAge_days=cumsum(days_since_prevObs),
maxAge_yrs=maxAge_days/365,
surv_01=as.numeric(doesSurvive),
logHoldfast=log(Holdfast)) %>%
filter(!is.na(doesSurvive))
surv.df %>% group_by(Management) %>%
summarise(nPlants=n_distinct(paste(Site, Patch, Plant)),
nMortalities=sum(!doesSurvive),
nSurvivals=sum(doesSurvive))
## # A tibble: 2 × 4
## Management nPlants nMortalities nSurvivals
## <fct> <int> <int> <int>
## 1 OA 72 59 24
## 2 TURF 53 17 57
ggplot(surv.df, aes(Holdfast, as.numeric(doesSurvive))) +
geom_jitter(shape=1, height=0.025, width=0.25) +
stat_smooth(method="glm", method.args=list(family="binomial")) +
scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20))
ggplot(surv.df, aes(Holdfast, as.numeric(doesSurvive), colour=Management)) +
geom_jitter(shape=1, height=0.025, width=0.25) +
stat_smooth(method="glm", method.args=list(family="binomial")) +
scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20))
# This is asking a lot of the data, but we can get reasonable looking curves
# across management types. The outcome isn't all that sensitive to priors
# at least within reason. Here's an option expecting greater variability among
# sites, then among plants, then among patches, and some nudged slopes.
library(brms)
out <- brm(bf(surv_01 ~ inv_logit(pSurvYr)^maxAge_yrs,
pSurvYr ~ logHoldfast * Management +
(1+logHoldfast * Management|Site/Patch/Plant),
nl=T),
data=surv.df, init=0, family="bernoulli", chains=4, cores=4,
prior=c(prior(normal(0, 1), nlpar="pSurvYr"),
prior(normal(-1, 1), nlpar="pSurvYr", coef="Intercept"),
prior(normal(1, 1), nlpar="pSurvYr", coef="ManagementTURF"),
prior(normal(0, 0.5), nlpar="pSurvYr", class="sd", lb=0,
group="Site"),
prior(normal(0, 0.1), nlpar="pSurvYr", class="sd", lb=0,
group="Site:Patch"),
prior(normal(0, 0.25), nlpar="pSurvYr", class="sd", lb=0,
group="Site:Patch:Plant")))
fitted.df <- expand_grid(logHoldfast=seq(-2.5,3.5,length.out=100),
Management=levels(surv.df$Management)) %>%
mutate(maxAge_yrs=1,
Holdfast=exp(logHoldfast),
Site=NA, Patch=NA, Plant=NA) %>%
bind_cols(posterior_epred(out, newdata=., nlpar="pSurvYr",
ndraws=4000, re_formula=NA) %>%
t %>% as_tibble) %>%
pivot_longer(starts_with("V"), names_to="iter", values_to="logit_p") %>%
group_by(logHoldfast, Holdfast, Management) %>%
summarise(logit_p_mn=mean(logit_p),
logit_p_lo=tidybayes::qi(logit_p, 0.9)[,1],
logit_p_hi=tidybayes::qi(logit_p, 0.9)[,2]) %>%
mutate(across(starts_with("logit_p"), inv_logit_scaled,
.names="{str_sub({.col}, 7, -1)}"))
fitted.df %>%
ggplot(aes(Holdfast, p_mn, colour=Management, fill=Management)) +
geom_ribbon(aes(ymin=p_lo, ymax=p_hi), alpha=0.5, colour=NA) +
geom_line() +
geom_jitter(data=surv.df, aes(y=surv_01),
shape=1, height=0.025, width=0.05) +
scale_x_continuous(trans="log", breaks=c(0.1, 1, 5, 20)) +
labs(x="Holdfast diameter (cm)", y="Annual survival probability")
fitted.df %>%
ggplot(aes(logHoldfast, logit_p_mn, colour=Management, fill=Management)) +
geom_ribbon(aes(ymin=logit_p_lo, ymax=logit_p_hi), alpha=0.5, colour=NA) +
geom_line() +
geom_jitter(data=surv.df, aes(y=logit_scaled(pmax(pmin(surv_01, 0.999), 0.001))),
shape=1, height=0.025, width=0.05) +
labs(x="log holdfast diameter (cm)", y="Annual survival probability")
We can link the holdfast size at t+1 to holdfast size at t x Management. Still to do: account for differences in number of days between monitoring surveys. Calculate seasonal or annual rates as possible.
grow.df <- Growth_Surv %>%
select(Site, Management, Patch, Monitoring,
Plant, Holdfast, Stipes, Total_Length) %>%
left_join(., survey.i) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
mutate(prevSurv=Monitoring-1) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
group_by(Site, Patch, Plant) %>%
mutate(doesSurvive=lead(Holdfast>0),
nextHoldfast=lead(Holdfast),
nextStipes=lead(Stipes),
nextLength=lead(Total_Length),
maxAge_days=cumsum(days_since_prevObs),
maxAge_yrs=maxAge_days/365) %>%
filter(doesSurvive) %>%
ungroup %>%
mutate(annHoldfastGrowth=(nextHoldfast-Holdfast)/days_to_nextObs*365,
annStipeGrowth=(nextStipes-Stipes)/days_to_nextObs*365,
annLengthGrowth=(nextLength-Total_Length)/days_to_nextObs*365)
ggplot(grow.df, aes(Holdfast, Holdfast+annHoldfastGrowth)) +
geom_abline(colour="grey") +
geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2)) +
scale_x_log10() + scale_y_log10()
ggplot(grow.df, aes(Holdfast, Holdfast+annHoldfastGrowth, colour=Management)) +
geom_abline(colour="grey") +
geom_point() +
stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
stat_smooth(method="lm", formula=y~x+I(x^2)) +
scale_x_log10() + scale_y_log10()
ggplot(grow.df, aes(Stipes, Stipes+annStipeGrowth)) +
geom_abline(colour="grey") +
geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2))
ggplot(grow.df, aes(Stipes, Stipes+annStipeGrowth, colour=Management)) +
geom_abline(colour="grey") +
geom_point() +
stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
stat_smooth(method="lm", formula=y~x+I(x^2))
ggplot(grow.df, aes(Total_Length, Total_Length+annLengthGrowth)) +
geom_abline(colour="grey") +
geom_point() + stat_smooth(method="lm", formula=y~x+I(x^2))
ggplot(grow.df, aes(Total_Length, Total_Length+annLengthGrowth, colour=Management)) +
geom_abline(colour="grey") +
geom_point() +
stat_smooth(aes(group=Site), method="lm", formula=y~x+I(x^2), se=F, size=0.25) +
stat_smooth(method="lm", formula=y~x+I(x^2))
We can estimate size distributions of recruits by management, and estimate the number of recruits as a function of the number or size of non-recruits in the plot. TODO: Account for the differences in time between monitoring periods to annualize somehow.
rcr.df <- Growth_Surv %>%
select(Site, Management, Patch, Monitoring,
Plant, Holdfast, Stipes, Total_Length) %>%
left_join(., survey.i) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
mutate(prevSurv=Monitoring-1) %>%
arrange(Site, Patch, Plant, Monitoring) %>%
group_by(Site, Patch, Plant) %>%
mutate(newRecruit=is.na(lag(Holdfast)),
maxAge_days=cumsum(days_since_prevObs),
maxAge_yrs=maxAge_days/365) %>%
group_by(Site, Patch, Monitoring) %>%
mutate(nPlants=sum(Holdfast>0),
totHoldfast=sum(Holdfast)) %>%
ungroup %>%
filter(newRecruit==TRUE) %>%
filter(maxAge_yrs < 1.5)
# Some bias in measurements: preference for 0.2, 10 cm
table(rcr.df$Holdfast)
##
## 0.1 0.2 0.6 0.8 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 8 60 1 1 8 6 6 5 3 6 1 3 1 10 1 6 1 1 3 3
## 18 20 22 23 24 28 30
## 2 1 2 1 1 1 2
ggplot(rcr.df, aes(Holdfast)) + geom_density() +
ggtitle("Recruit size distr.") +
scale_x_continuous(trans="log1p")
ggplot(rcr.df, aes(Holdfast, colour=Management)) +
geom_density(aes(group=Site), size=0.25) +
geom_density(size=1) +
ggtitle("Recruit size distr.") +
scale_x_continuous(trans="log1p")
ggplot(rcr.df, aes(maxAge_yrs)) + geom_histogram() +
ggtitle("Recruit max possible age")
ggplot(rcr.df, aes(Holdfast)) + geom_histogram() +
scale_x_continuous(trans="log1p") + facet_grid(Management~Monitoring)
ggplot(rcr.df, aes(factor(Monitoring), fill=Management)) +
geom_bar(position="dodge") +
ggtitle("Number of recruits decreases with time")
ggplot(rcr.df, aes(Management, fill=factor(Monitoring))) +
geom_bar(position="fill") +
ggtitle("Number of recruits decreases with time, maybe more strongly in TURF")
rcr.df %>% group_by(Site, Management, Patch, Monitoring, nPlants) %>%
summarise(nRecruits=sum(newRecruit)) %>%
ggplot(aes(nPlants-nRecruits, nRecruits)) +
geom_jitter(width=0.2, height=0.2) +
stat_smooth(method="lm", formula=y~log1p(x)) +
scale_y_continuous(trans="log1p") +
xlab("Number of non-recruits")
rcr.df %>% group_by(Site, Management, Patch, Monitoring, nPlants) %>%
summarise(nRecruits=sum(newRecruit)) %>%
ggplot(aes(nPlants-nRecruits, nRecruits, colour=Management)) +
geom_jitter(width=0.2, height=0.2) +
stat_smooth(aes(group=Site), method="lm", formula=y~log1p(x), se=F, size=0.25) +
stat_smooth(method="lm", formula=y~log1p(x), size=1) +
scale_y_continuous(trans="log1p") +
xlab("Number of non-recruits")
rcr.df %>% group_by(Site, Management, Patch, Monitoring, totHoldfast) %>%
summarise(nRecruits=sum(newRecruit),
HoldfastRcr=sum(newRecruit*Holdfast)) %>%
ggplot(aes(totHoldfast-HoldfastRcr, nRecruits)) +
geom_jitter(width=0.2, height=0.2) +
stat_smooth(method="lm", formula=y~log1p(x)) +
scale_y_continuous(trans="log1p") +
xlab("Total holdfast diameter of non-recruits")
rcr.df %>% group_by(Site, Management, Patch, Monitoring, totHoldfast) %>%
summarise(nRecruits=sum(newRecruit),
HoldfastRcr=sum(newRecruit*Holdfast)) %>%
ggplot(aes(totHoldfast-HoldfastRcr, nRecruits, colour=Management)) + geom_jitter() +
stat_smooth(aes(group=Site), method="lm", formula=y~log1p(x), se=F, size=0.25) +
stat_smooth(method="lm", formula=y~log1p(x), size=1) +
scale_y_continuous(trans="log1p") +
xlab("Total holdfast diameter of non-recruits")
# The recruit size distribution seems to change with time-since-harvest.
# Mainly small recruits at the start, then larger recruits at TURF sites as
# populations become established... These should only be *new* recruits since
# the last survey. Differences in thinning / establishment by management? OA
# tend to have more and smaller recruits after the initial pulse.
# How does this interact with higher survival rates and faster growth rates in TURF?
# Lower establishment, but better survival and growth once established?
ggplot(rcr.df, aes(days_since_harvest, Holdfast)) +
geom_jitter(width=10, height=0.1) +
stat_smooth(span=1.5) +
ggtitle("Recruit size vs. time since clearance") +
scale_y_continuous(trans="log1p")
ggplot(rcr.df, aes(days_since_harvest, Holdfast, colour=Management)) +
geom_jitter(width=10, height=0.1) +
stat_smooth(aes(group=Site), span=1.5, se=F, size=0.25) +
stat_smooth(span=1.5) +
ggtitle("Recruit size vs. time since clearance") +
scale_y_continuous(trans="log1p")
ggplot(rcr.df, aes(Holdfast, colour=factor(Monitoring))) +
geom_density(aes(group=paste(Site, Monitoring)), adjust=1.5) +
ggtitle("Recruit size vs. time since clearance") +
scale_y_continuous(trans="log1p") +
facet_wrap(~Management)
rcr.df %>% group_by(Site, Management, Patch, Monitoring, days_since_harvest) %>%
summarise(nRecruits=sum(newRecruit)) %>%
ggplot(aes(days_since_harvest, nRecruits)) +
geom_jitter(width=10, height=0.1) +
stat_smooth(span=1.5) +
ggtitle("Number of recruits vs. time since clearance") +
scale_y_continuous(trans="log1p")
rcr.df %>% group_by(Site, Management, Patch, Monitoring, days_since_harvest) %>%
summarise(nRecruits=sum(newRecruit)) %>%
ggplot(aes(days_since_harvest, nRecruits, colour=Management)) +
geom_jitter(width=10, height=0.1) +
stat_smooth(aes(group=Site), span=1.5, se=F, size=0.25) +
stat_smooth(span=1.5) +
ggtitle("Number of recruits vs. time since clearance") +
scale_y_continuous(trans="log1p")
We can fit a nice curve predicting the probability of reproductive blades given holdfast diameter, but do we have any data linking reproductive blades to recruits? That is, what we need is \(F(z',z) \sim p_{repro}(z) f_{spores}(z) p_{estab} f_{rcrSize}(z')\). We have an estimate for \(p_{repro}(z)\), which is the probability of having reproductive blades given a holdfast size \(z\). We have part of \(f_{rcrSize}(z')\), which is the expected recruit size distribution at time \(t+1\). What we’re missing is how many spores might be produced by each reproductive plant, and how likely those spores are to establish.
We could potentially lump several processes and calculate directly the number of recruits from the data above based on a presumed maximum established recruits, with the actual as a function of \(N\).
morph.df %>%
filter(!is.na(Reproductive)) %>%
droplevels() %>%
ggplot(aes(Holdfast, Reproductive)) +
geom_jitter(shape=1, height=0.025, width=0.5) +
stat_smooth(method="glm", method.args=list(family="binomial")) +
scale_x_continuous(trans="log1p")
morph.df %>%
filter(!is.na(Reproductive)) %>%
droplevels() %>%
ggplot(aes(Holdfast, Reproductive, colour=Management)) +
geom_jitter(shape=1, height=0.025, width=0.5) +
stat_smooth(method="glm", method.args=list(family="binomial")) +
scale_x_continuous(trans="log1p")